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We numerically study spin turbulence in a two-dimensional trapped spin-1 spinor Bose-Einstein 
condensate, focusing on the energy spectrum. The spin turbulence in the trapped system is generated 
by instability of the helical structure of the spin density vector in the initial state. Our numerical 
calculation finds that in the trapped system the spectrum of the spin-dependent interaction energy 
in the ferromagnetic case exhibits a —7/3 power law, which was confirmed in a uniform system by 
our previous study. The relation between the —7/3 power law and the motion of the spin density 
vector is discussed by investigating the orbits of dynamical variables in the spin space. 
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I. INTRODUCTION 

Turbulence in classical fluids has been studied for a 
long time. It is an important theme in modern physics, 
but is not yet understood sufficiently [1]. Great progress 
in the study of turbulence was made by Kolmogorov [2] , 
who suggested that the spectrum of the kinetic energy 
exhibits a —5/3 power law in fully developed isotropic 
homogeneous turbulence. This is called the Kolmogorov 
—5/3 power law. 

Recently, quantum turbulence (QT) in one-component 
atomic Bose-Einstein condensates (BECs) has been in- 
vestigated. Some theoretical and numerical studies show 
that the Kolmogorov —5/3 power law appears in QT [3]. 
QT in an atomic BEC has been experimentally realized 
by Henn etal. [4], but the Kolmogorov —5/3 power law 
has not yet been observed. 

Multi-component atomic BECs can be realized [5], 
described by more than two macroscopic wave func- 
tions, which means that the degrees of freedom in 
multi-component BECs is larger than in one-component 
BECs. Thus, various behaviors that do not appear in 
one-component BEC are expected to occur in this sys- 
tem. The hydrodynamical instability in two-component 
BECs, such as Rayleigh-Taylor, Kelvin-Helmholtz and 
Richtmyer-Meshkov instability, has been studied using 
the Gross-Pitaevskii (GP) equation [6-8]. The turbu- 
lence in a two-component BEC has been numerically and 
theoretically investigated, and a tangle of two kinds of 
quantized vortices was found to form in the system [9, 10] . 

Another multi-component BEC system is the spinor 
BEC, comprised of condensates corresponding to differ- 
ent hyperfine states of an atom. We investigated the 
turbulence in a spinor BEC, namely spin turbulence, in 
our previous work [11]. In the prior research, we gen- 
erated the spin turbulence by a counterflow instability, 
finding that, in the ferromagnetic case, the spectrum of 
the spin-dependent interaction energy exhibits a —7/3 
power law. This power law can be understood by scaling 
analysis with some assumptions. One assumption is that 
the total density of the condensate is uniform. However, 



in the actual experiments, the condensates are trapped 
by a harmonic potential that is not uniform. Therefore, 
whether the —7/3 power law appears in the trapped sys- 
tem is an important issue. 

In this paper, we report the numerical observation of 
the —7/3 power law in a trapped spin-1 spinor BEC. Our 
study is based on a numerical calculation of the GP equa- 
tion. The spin turbulence in the trapped system arises 
from the instability of the helical structure of the spin 
density vector in the initial state. This helical structure 
has been experimentally realized by Vengalattore etal. 
[12]. We find that, in the ferromagnetic case, the —7/3 
power law also appears in the trapped system if the sys- 
tem size is much larger than the spin coherence length. 
This paper is organized as follows. Section II describes 
the formulation. The results of the numerical simulation 
of the GP equation are described in Sec. III. Here we 
show how the helical structure of the spin density vector 
becomes unstable to the development of spin turbulence 
and the —7/3 power law appears in the spin-dependent 
interaction energy. In Sec. IV, in order to understand 
the relation between the energy spectra and the motion 
of the spin density vector, we discuss the characteristic 
time of the spin, the orbits in the spin space and the size 
of the BEC cloud for sustaining the power law. Section 
V is devoted to conclusions. 



II. FORMULATION 

A. Gross— Pitaevskii equation 

We consider a trapped spin-1 spinor BEC at zero tem- 
perature, which is well described by the macroscopic 
wave functions ijj m (m = 1, 0, —1). Then, the wave 
functions ij) m obey the GP equation [13, 14]: 
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where M and V are the mass of a particle and the trap- 
ping potential. In our study, the potential is given by 
V = Muj 2 (x 2 + y 2 )/2 + Muj 2 z 2 /2. The total density n 
and the spin density vector Si (i = x,y,z ) are given 

by n = EL-i l^m| 2 and Si = El,n=-lfc( S i)^n, 
where (Si) mn are the spin-1 matrices. The parameters 
Co and c\ are the coefficients of the spin-independent and 
dependent interactions, which are expressed by 47rft 2 (ao + 
2<22 )/3M and 47rh 2 (a2 — ao)/3M with s wave scattering 
lengths ao and g&2. 

The total energy E is given by 

m =-l 

+^ y n 2 dr + c ± j s 2 dr . (2) 

The spin-dependent interaction energy is the last term 
with the coefficient c\ of the right hand side of Eq. (2), 
which is characteristic of the spinor BEC. The ground 
state in a uniform system without a magnetic field is 
ferromagnetic for c\ < and polar for c\ > 0. This 
interaction exchanges the particles between the different 
condensates, so the particle number in each component 
is not conserved. Thus, the instability and turbulence in 
spinor BECs are expected to have properties that one- 
and two-component BECs do not. 

We consider that the condensate is strongly confined 
in the z direction (u <C uo z ). In this case, we can approx- 
imately separate the degrees of freedom of the macro- 
scopic wavefunctions as ijj m (x,y,z,t) = ip rn (x^y^t)f(z). 
We assume f(z) = (1/2™^ ^e^ 2 / 4 ^* with a hz = 
(H/2Muj z ) 1 ^ 2 . Then, the GP equation (1) is reduced to 
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ih—lpm = (-^V 2 + V)VVi+C nVVn+Ci ^ S'Smn^n 
^ n=— 1 

(3) 

with V 2 = d 2 /dx 2 +d 2 /dy 2 , V = Muj 2 (x 2 + y 2 )/2, c = 
c /2y/7ra hzi ci = c\j2^a hz , n = Y^m=-i \^™\ 2 and 
= El,n=-iCK)mni- We consider experiments 
with 87 Rb, which exhibits a ferromagnetic interaction. 
Thus, we use M = 1.42 x 10" 25 kg, a = 5.39 x 10" 9 m, 
a 2 = 5.31 x 10" 9 m, TV = 3 x 10 5 , u = 2tt x 20 /s and 
uj z = 2tt x 600 /s. 

B. Numerical method 

We use the Crank-Nicholson method to numerically 
calculate the GP equation (3). The coordinate is nor- 
malized by the length = (h/2Muj) 1 ^ 2 ~ 1.72/im and 
the box size is 80*2^ x 80a/!. Space in the x and y direc- 
tions is discretized into 1024 x 1024 bins. The time is 
normalized by the frequency uo of the trapping potential 
in the x and y directions. 




FIG. 1: (Color online) Profiles of each component of the spin 
density vector in the initial state. The box size is 40a^ x 40a^, 
where a h is (h/2Muj) 1/2 - 1.72/im. 



C. Initial state 

We use an initial state with a helical structure of the 
spin density vector to obtain the spin turbulence in the 
trapped spin-1 spinor BEC. In the following, we show the 
mathematical expression of this state. The macroscopic 
wave functions i\) = ((/>, 0, 0) has spin density vector s = 
|0| 2 e 2 , where ej(j = x,y,z) is the unit vector in the 
j direction. Here, can be obtained by the imaginary 
time step of Eq. (3). By multiplying this wave function 
by the rotation matrix U = e - %OLS z e -^S Ve -%-iS z m ^ 
spin space, we obtain 

/e— cos 2 f\ 
Uo = fc-* Tisin/? , (4) 
V?-i/ V e-sin 2 f / 

where a, /3, 7 are the Euler angles. Then, the spin 
density vector is expressed by s = |(/>| 2 (sin/3 cosa e x + 
sin/3 sina e y + cos/3 e z ). Therefore, using a = 7r/2, 
f3 = khX and 7 = 0, the spin density vector has a he- 
lical structure with wave number kh- In this paper, we 
consider the case of kh ~ 1.22 //am. Thus, Eq. (4) be- 
comes 

^1 \ _ /-icos 2 (k h x/2)\ 
#0 = -^sin(k h x) . (5) 
^-1/ \ isin 2 (k h x/2) J 

Figure 1 shows the profile of each component of the spin 
density vector in the initial state of Eq.(5). The heli- 
cal structure has been experimentally realized by Ven- 
galattore etal. [12], where they prepared the structure by 
means of a magnetic field, investigating how the structure 
becomes unstable to changes into some disordered state 
through observing the spin density vector. In our calcu- 
lations, we add some small white noise to the initial state 
of Eq. (5). The noise causes the particle number of each 
component to fluctuate by 0.1-0.3%. This is consistent 
with experimental results [15]. 
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(a) uot = 



(b) out = 5 





(c) cot 60 



(d) out = 180 





FIG. 2: (Color online) Profiles of the spin density vector at 
uot — 0, 5, 60, 180. The shading of the arrows denotes the 
amplitude of the spin density vector. The box size is 40a^, x 
40a^. 



D. Spectrum of spin-dependent interaction energy 

We derive an expression for the spectrum of the spin- 
dependent interaction energy. The spin-dependent inter- 
action energy E s is given by 



E, 



dr 
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with f = (x,y). We expand the spin density vector s(f) 
with plane waves:s(f) = s(fe)e* fc *' r with k — (k xi k y ). 
Then, the spin- dependent interaction energy E s is repre- 
sented by s(k) as 



E s = c -^Y,\ s Ck)\ 2 , 

k 
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where A is the area of the system. Therefore, the en- 
ergy spectrum of the spin-dependent interaction energy 
is given by 



E.(k) 
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/c<|fei|<fe+A/c 

where Ak is 2tt/L for a system size L. 

III. NUMERICAL RESULTS 

We show our numerical results with the GP equation 
(3) starting from the initial state of Eq. (5) [16]. Our 



calculation finds that the helical structure in the trapped 
system leads to spin turbulence and the spectrum of 
the spin-dependent interaction energy exhibits the —7/3 
power law. 

The helical structure of the spin density vector in the 
trapped system is unstable, generating the spin turbu- 
lence. Figure 2 (a) shows the profile of the spin density 
vector in the initial state. This helical structure immedi- 
ately collapses as shown in Fig. 2 (b). As time goes by, 
the spin density vector is significantly disturbed, point- 
ing in various directions, as shown in Figs. 2 (c) and (d). 
Thus, spin turbulence in the trapped system is realized. 

The spectrum of the spin-dependent interaction energy 
exhibits the —7/3 power law in the spin turbulence in the 
trapped system. Figures 3 show the time development 
of the spectrum, which is numerically calculated by Eq. 
(8). In the initial state, the spectrum in Fig. 3 (a) has a 
sharp peak at kan ~ 2 corresponding to the wave number 
kh of the helical structure. As the instability develops, 
modulation of the spin density vector with various wave 
numbers is excited. Thus, as shown in Fig. 3 (b), the 
spectrum loses the sharp peak, letting the energy flow 
into both the low and high wave number regions. As 
time passes, the spin density vector is significantly dis- 
turbed and the system exhibits spin turbulence. Then, 
the spectrum develops to exhibit the —7/3 power law in 
the wave number region krF < k < /c s , as shown in Fig. 
3 (c). Here &tf an d k s are the wave numbers corre- 
sponding to the Thomas-Fermi radius Rtf and the spin 
coherence length £ s = h/ ^/2M\c\ \uq. This —7/3 power 
law in the spectrum continues at least until ut = 180. 
Figure 3 (d) shows the spectrum at ut = 180. 



A. 



IV. DISCUSSION 

Motion of the spin density vector and the 
power law 



-7/3 



The appearance of the —7/3 power law is accompanied 
by characteristic behavior of the spin density vector. We 
should note the time scale of the motion as well as the 
randomness of the spin density vector. After the ini- 
tial helical structure becomes unstable, the spin vectors 
exhibit relatively large-scale, slow motion. The helical 
structure then disappears and the vectors become disor- 
dered. As the —7/3 power law appears, the vectors begin 
to perform rapid, fine oscillations around random direc- 
tions; the motion of each spin vector appears to be frozen 
to a random direction. 

This can be understood qualitatively by the follow- 
ing discussion. In our previous study [11], we performed 
Kolmogorov-type scaling analysis for the equation of mo- 
tion of the spin density vector and thus obtained the 
—7/3 power law. In the wave number region sustaining 
the —7/3 power law, the characteristic time t s (k) of the 
mode of wave number k is expressed by 

t s (k) -e- 1 / 3 ^ 3 / 4 , (9) 
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FIG. 3: (Color online) Spectrum of the spin-dependent interaction energy at out — 0, 5, 60, 180. k s and krF are the wave 
numbers corresponding to the spin coherence length £ s — h/^/2M\a \no and the Thomas-Fermi radius Rtf, where no is the 
total density at the center in the initial state. Their expressions are k s = 27r/£ s and krF — 2tt/Rtf- 



where e is the energy flux flowing from low to high wave 
numbers. Equation (9) is obtained from e ~ k~Hj 3 , 
which means that the energy flux is independent of the 
wave number [11]. Equation (9) shows that the charac- 
teristic time t s shortens when the wave number is high. 

As the spin density vectors become disordered and the 
high k modes (small scale) are excited, the characteristic 
time scale becomes short, leading to the rapid, fine oscil- 
lations. On the other hand, the motion with low k modes 
(large scale) has a long characteristic time t s . Thus, the 
rough direction of the spin density vector at each position 
is largely fixed and changes only slowly. 

Thus the motion of the spin density vector is closely 
related to the advent of the —7/3 power law. 



B. Orbits in spin space 

In order to reveal these properties, it is useful to 
consider the orbits of the dynamical variables in phase 
space, which is an approach often used in the re- 
search field of nonlinear dynamics. Figure 4 shows 
the orbits of the normalized spin density vectors in 
the (s x a\/N, s y a\/N, s z a\/N) space for different time 
stages. The two orbits correspond to the spin density vec- 
tor at the positions (0, 0) and (0, Rtf/2) inside the BEC 
cloud. The amplitude of the spin density vector depends 
on the position inside the BEC cloud and is maximum 



at the center; such a dependence of the amplitude is re- 
flected in the orbits of the spin at different positions. We 
discuss the time development of the orbits of Fig. 4 with 
reference to Fig. 2 and Fig. 3. The periodic orbits in the 
early stage (Fig. 4 (a)) indicate some coherent motion of 
the initial helical structure. Some of the orbits become 
unstable (Fig. 4 (b)) after the helical structure become 
unstable. Then, the spin density vectors start to become 
disordered (Fig. 2 (b)) and the energy spectrum becomes 
broad but does not yet exhibit the —7/3 power law (Fig. 
3(b)). After a time, the two orbits become separated, as 
shown in Fig. 4(c), and independently localized in some 
places in the spin space, generating rapid, fine oscilla- 
tion in those areas. During this process the spin density 
vectors become disordered and frozen (Fig. 2(c)) and 
the energy spectrum exhibits the —7/3 power law (Fig. 
3(c)). This behavior continues until cut = 180, while the 
orbits appear to shrink from Fig. 4 (b) through (c) to 
(d). After the power law is established, the orbits move 
very slowly in the spin space, as shown from Fig. 4 (c) 
to (d). 

In order to measure the shrinkage quantitatively, we 
introduce the size I of the orbit in the spin space, where 
the size I at time r = uot is defined by the maximum 
distance between two points in an orbit during r tor + 2. 
Figure 5 shows the sizes l(r) of orbits at three different 
positions (0,0), (0, Rtf/2) and (0,3.Rtf/4) inside the 
BEC cloud. The sizes of the orbits are found to fluctuate 
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FIG. 4: (Color online) Orbits of the spin density vector at two positions (0,0) (red filled squares) and (0,Rtf/2) (blue empty 
squares) inside the BEC cloud for the periods (a) ut = 0-2, (b) 5-7, (c) 60-62 and (d) 180-182. The spin density vector is 
normalized by N/a\. 



and shrink as time passes; the size is generally smaller at 
outer positions in the cloud, presumably because of the 
smaller amplitude of the spin density vector. 

It is interesting that the freezing of the separated spin 
orbits can determine the power law of the energy spec- 
trum. 



C. Size of condensate and the —7/3 power law 

We discuss the relation between the size of the con- 
densate and the —7/3 power law. We confirm that if 
the Thomas-Fermi radius Rtf is not sufficiently larger 
than the spin coherence length £ s , the spectrum of the 
spin-dependent interaction energy does not exhibit the 
—7/3 power law. This is confirmed by a numerical cal- 
culation with Rtf /is ~ 5. For the calculations in this 
paper, the condition Rtf Us ^ 10 is satisfied. Therefore, 
in the experiments, the condition Rtf/ is > 10 could be 
necessary for observing the —7/3 power law. 



The spin turbulence in the trapped system is generated 
by instability of the helical structure of the spin density 
vector. Our numerical calculation finds that the spec- 
trum of the spin-dependent interaction energy exhibits 
the —7/3 power law, which is consistent with our pre- 
vious result for a uniform system [11]. We discuss the 
relation between the —7/3 power law and the motion of 
the spin density vector. 

We believe that the —7/3 power law can be experi- 
mentally observed because of the following reasons. As 
discussed in Sec. IV C, the spin turbulence in a large 
trapped system (Rtf /is ~ 10) sustains this power law 
for a long time. All the components of the spin density 
vector have been experimentally observed by a phase con- 
trast imaging method [12, 15]. The expressions for the 
spectrum of the spin-dependent interaction energy show 
that it is possible to obtain a spectrum if the spin density 
vector is observed everywhere. This is because Eq. (8) 
contains only the Fourier component of the vector s(k). 
Therefore, it is possible to observe this power law. 



CONCLUSIONS 
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FIG. 5: (Color online) Time development of the size of the 
orbits of the spin density vector in the spin space of Fig. 
4. The three curves correspond to the three positions (0,0), 
(0, Rtf/2) and (0, 3Rtf/4) inside the BEC cloud. The size / 
at time r = out refers to the maximum distance between two 
points in an orbit during r to r + 2. 
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